
#### Preliminaries ####

R.version$version.string

# clear workspace
rm(list=ls())

# load packages
library(here)
library(tidyverse)
library(rworldmap) # maps
library(magrittr)
library(reshape2)
library(ggpubr)
library(leaflet)
library(plotly)
library(gapminder)
library(countrycode)
library(ggrepel)
library(ggthemes)
library(modelsummary)


ggplot_theme <- theme(axis.text=element_text(size=16),
                      axis.title=element_text(size=16,face="bold"))

#---------------------------------------------------#
#--------------------load data ---------------------#
#---------------------------------------------------#

# 1. Load V-Dem data (version 12)

# set you working directory here #

#setwd()

vdem_full_v13 <- readRDS("data/vdem/V-Dem-CY-Full+Others-v13.rds")


vdem_subset <- vdem_full_v13 %>%
  dplyr::select(country_name, year, country_id, starts_with("v2x_polyarchy"),
                starts_with("v2cafres"), starts_with("v2cafexch"), 
                starts_with("v2cainsaut"), starts_with("v2casurv"), 
                starts_with("v2clacfree"))

vdem <- vdem_full_v13

#---------------------------------------------------#
#---- Build Function for getting Episodes of X -----#
#---------------------------------------------------#

# load the relevant functions into the environment

source("R/get_episode_function.R")
source("R/get_episode_wo_CI_function.R")
source("R/plot_all_episode.R")
source("R/plot_episodes_function.R")


#---------------------------------------------------#
#---------Episodes of Academic Freedom Increase-----#
#---------------------------------------------------#

episode_data <- get_episode(data = vdem)
episode_data_wo_CI <- get_episode_wo_CI(data = vdem)


################################################################################
########################## FIGURE E2 Appendix ##################################
################################################################################

#### Plot AFI decreases ####

episode_data <- episode_data %>%
  group_by(country_name) %>%
  mutate(lag_v2xca_academ = dplyr::lag(v2xca_academ)) %>%
  ungroup() %>%
  dplyr::select(country_id, year, starts_with("increase_"), starts_with("decline"), lag_v2xca_academ)

vdem_full_v13 <- vdem_full_v13 %>%
  left_join(episode_data, by = c("country_id", "year"))

vdem_subset <- vdem_full_v13 %>%
  dplyr::select(country_name, country_id, year, country_text_id, starts_with("v2xca_academ"), 
                starts_with("increase_"), starts_with("decline"), lag_v2xca_academ, v2cauni, starts_with("gap"), v2x_regime, 
                codingstart, codingend) %>%
  filter(year >=1900) %>%
  
  # prepare dataset for plotting by constructing a full.df with gaps and NAs for the gaps. 
  
  dplyr::arrange(country_text_id, year) %>%
  dplyr::group_by(country_id) %>%
  # make codingstart 1900 or first year thereafter
  dplyr::mutate(codingstart2 = min(hablar::s(ifelse(!is.na(v2x_regime), year, NA))),
                # tag original sample for later use
                origsample = 1) %>%
  # we need to deal with gaps in v-dem coding
  # this balances the dataset
  plm::make.pconsecutive(balanced = TRUE, index = c("country_id", "year")) %>%
  dplyr::group_by(country_id) %>%
  # this fills missing variables we need that are constant within countries
  tidyr::fill(c(country_text_id, country_name, codingend, gapstart1, gapend1, gapstart2, gapend2,
                gapstart3, gapend3)) %>%
  tidyr::fill(c(country_text_id, country_name,codingend, gapstart1, gapend1, gapstart2, gapend2,
                gapstart3, gapend3), .direction = "up")  %>%
  # here we need to recode the gaps as only during the period prior to and during the gap (for our "uncertain" variables)
  dplyr::mutate(gapstart = ifelse(year <= gapend1, gapstart1, NA),
                gapend = ifelse(year <= gapend1, gapend1, NA),
                gapstart = ifelse(!is.na(gapend2) & year > gapend1 & year <= gapend2, gapstart2, gapstart),
                gapend = ifelse(!is.na(gapend2) & year > gapend1 & year <= gapend2, gapend2, gapend),
                gapstart = ifelse(!is.na(gapend3) & year > gapend2 & year <= gapend3, gapstart3, gapstart),
                gapend = ifelse(!is.na(gapend3) & year > gapend2 & year <= gapend3, gapend3, gapend))

# plotting data #

f1 <- ggplot() +
  geom_line(subset(vdem_subset),
            mapping = aes(x = year, y = v2xca_academ, group = country_text_id, color = as.factor(decline_episode)), size = 0.2, alpha = 0.2) +
  geom_line(subset(vdem_subset, decline_episode==1), 
            mapping = aes(x = year, y = v2xca_academ, group = decline_episode_id, color = as.factor(decline_episode)), size = 0.5) +
  theme_pubr() + 
  scale_color_manual(values = c("#d3d3d3", "#920050", "#d3d3d3"), labels = c("No Episode", "Episode", "NA")) +
  scale_x_continuous(breaks = seq(round(min(vdem_subset$year) / 10) * 10, round(max(vdem_subset$year) / 10) * 10, 10)) +
  labs(y = "Academic Freedom Index",
       color = "", x = "Year",
       size = "",
       shape = "",
       title = "",
       subtitle = "",
       caption = "") +
  ggplot_theme

f2 <- ggplot() +
  geom_line(subset(vdem_subset),
            mapping = aes(x = year, y = v2xca_academ, group = country_text_id, color = as.factor(increase_episode)), size = 0.2, alpha = 0.2) +
  geom_line(subset(vdem_subset, increase_episode==1), 
            mapping = aes(x = year, y = v2xca_academ, group = increase_episode_id, color = as.factor(increase_episode)), size = 0.5) +
  theme_pubr() + 
  scale_color_manual(values = c("#d3d3d3", "#FFB81C", "#d3d3d3"), labels = c("No Episode", "Episode", "NA")) +
  scale_x_continuous(breaks = seq(round(min(vdem_subset$year) / 10) * 10, round(max(vdem_subset$year) / 10) * 10, 10)) +
  labs(y = "Academic Freedom Index",
       color = "", x = "Year",
       size = "",
       shape = "",
       title = "",
       subtitle = "",
       caption = "") +
  ggplot_theme

vdem_subset_sum <- vdem_subset %>%
  group_by(year) %>%
  summarize(number_of_countries = sum(v2cauni, na.rm = T), 
            number_of_decliners = sum(decline_episode, na.rm = T), 
            number_of_increasers = sum(increase_episode, na.rm = T)) %>%
  mutate(number_of_countries = number_of_countries - number_of_decliners - number_of_increasers) %>%
  pivot_longer(!year, names_to = "type_country_nr", values_to = "count")

f3 <- ggplot(vdem_subset_sum, aes(x = year, y= count, color = type_country_nr, fill = type_country_nr)) +
  geom_bar(stat="identity", position="stack", width = 0.7) +
  theme_pubr() + 
  scale_color_manual(values = c("#d3d3d3", "#920050", "#FFB81C"),  labels = c("No. of Countries", "No. of Decliners", "No. of Increasers")) +
  scale_fill_manual(values = c("#d3d3d3", "#920050", "#FFB81C"), labels = c("No. of Countries", "No. of Decliners", "No. of Increasers")) +
  scale_x_continuous(breaks = seq(round(min(vdem_subset_sum$year) / 10) * 10, round(max(vdem_subset_sum$year) / 10) * 10, 10)) +
  labs(y = "Countries in Episode",
       color = "", x = "Year",
       size = "",
       fill = "",
       shape = "",
       title = "",
       subtitle = "",
       caption = "") +
  ggplot_theme

decline_increase_overview <- cowplot::plot_grid(f1+  theme(legend.position = "none"),f2 +  theme(legend.position = "none"), f3,
                                                ncol=1,
                                                align = "v",
                                                rel_heights = c(1.1,1.1,1))

#---------------------------------------------------#

ggsave(here("figures/decline_increase_overview_w_CI.pdf"), plot = decline_increase_overview, height = 20, width = 20)
ggsave(here("figures/decline_increase_overview_w_CI.png"), plot = decline_increase_overview, height = 20, width = 20)


################################################################################
########################## FIGURE E3 Appendix ##################################
################################################################################

#### Waves Argument ####

vdem_subset_waves <- vdem_subset_sum %>%
  pivot_wider(id_cols = year, names_from = type_country_nr, values_from = count) %>%
  dplyr::mutate(net_change = number_of_increasers - number_of_decliners, 
                roll_change = zoo::rollmean(net_change, 5, align = "center", na.pad = TRUE))

ggplot(vdem_subset_waves, aes(x = year, y= net_change)) +
  geom_bar(stat="identity", position="stack", width = 0.7) +
  geom_line(aes(y = roll_change), color = "#920050", size = 1.05) +
  theme_pubr() + 
  scale_x_continuous(breaks = seq(round(min(vdem_subset_sum$year) / 10) * 10, round(max(vdem_subset_sum$year) / 10) * 10, 10)) +
  labs(y = "Number of Growth Episodes -\nNumber of Decline Episodes",
       color = "", x = "Year",
       size = "",
       fill = "",
       shape = "",
       title = "",
       subtitle = "",
       caption = "") +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=16,face="bold"))

ggsave(here("figures/net_change_w_CI.pdf"),height = 10, width = 20, units = "cm", dpi = 800)
ggsave(here("figures/net_change_w_CI.png"),height = 10, width = 20, units = "cm", dpi = 800)


################################################################################
########################## FIGURE E4 Appendix ##################################
################################################################################

p1 <- plot_all(data.df = episode_data)
p2 <- ERT::plot_all()


cowplot::plot_grid(p1, p2, 
                   ncol=1,
                   align = "v",
                   rel_heights = c(1.1,1.1), 
                   labels = c('Academic Freedom', 'Political Regimes'))

ggsave("figures/Overview_Episodes_with_ERT_data.png", height =25, width = 30, units = "cm", dpi = 800)
